home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus 2000 #4 / Amiga Plus CD - 2000 - No. 4.iso / PowerPC / Dev / PPCRelease / Libmfd / e_remainder.c < prev    next >
Encoding:
C/C++ Source or Header  |  1997-01-02  |  1.7 KB  |  78 lines
  1.  
  2. /* @(#)e_remainder.c 1.3 95/01/18 */
  3. /*
  4.  * ====================================================
  5.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  6.  *
  7.  * Developed at SunSoft, a Sun Microsystems, Inc. business.
  8.  * Permission to use, copy, modify, and distribute this
  9.  * software is freely granted, provided that this notice 
  10.  * is preserved.
  11.  * ====================================================
  12.  */
  13.  
  14. /* __ieee754_remainder(x,p)
  15.  * Return :                  
  16.  *     returns  x REM p  =  x - [x/p]*p as if in infinite 
  17.  *     precise arithmetic, where [x/p] is the (infinite bit) 
  18.  *    integer nearest x/p (in half way case choose the even one).
  19.  * Method : 
  20.  *    Based on fmod() return x-[x/p]chopped*p exactlp.
  21.  */
  22.  
  23. #include "fdlibm.h"
  24.  
  25. #ifdef __STDC__
  26. static const double zero = 0.0;
  27. #else
  28. static double zero = 0.0;
  29. #endif
  30.  
  31.  
  32. #ifdef __STDC__
  33.     double __ieee754_remainder(double x, double p)
  34. #else
  35.     double __ieee754_remainder(x,p)
  36.     double x,p;
  37. #endif
  38. {
  39.     int hx,hp;
  40.     unsigned sx,lx,lp;
  41.     double p_half;
  42.  
  43.     hx = __HI(x);        /* high word of x */
  44.     lx = __LO(x);        /* low  word of x */
  45.     hp = __HI(p);        /* high word of p */
  46.     lp = __LO(p);        /* low  word of p */
  47.     sx = hx&0x80000000;
  48.     hp &= 0x7fffffff;
  49.     hx &= 0x7fffffff;
  50.  
  51.     /* purge off exception values */
  52.     if((hp|lp)==0) return (x*p)/(x*p);     /* p = 0 */
  53.     if((hx>=0x7ff00000)||            /* x not finite */
  54.       ((hp>=0x7ff00000)&&            /* p is NaN */
  55.       (((hp-0x7ff00000)|lp)!=0)))
  56.         return (x*p)/(x*p);
  57.  
  58.  
  59.     if (hp<=0x7fdfffff) x = __ieee754_fmod(x,p+p);    /* now x < 2p */
  60.     if (((hx-hp)|(lx-lp))==0) return zero*x;
  61.     x  = fabs(x);
  62.     p  = fabs(p);
  63.     if (hp<0x00200000) {
  64.         if(x+x>p) {
  65.         x-=p;
  66.         if(x+x>=p) x -= p;
  67.         }
  68.     } else {
  69.         p_half = 0.5*p;
  70.         if(x>p_half) {
  71.         x-=p;
  72.         if(x>=p_half) x -= p;
  73.         }
  74.     }
  75.     __HI(x) ^= sx;
  76.     return x;
  77. }
  78.